package com.bravo.other.geoid;

import java.util.ArrayList;
import java.util.List;

/**
 * Advanced Geoid Height Calculator
 * Contains complete spherical harmonic calculation and EGM2008 model implementation
 * 
 * @author bravo
 */
public class AdvancedGeoidCalculator {
    
    // Main spherical harmonic coefficients for EGM2008 model
    private static final double[][] C_COEFFICIENTS = {
        // n=0
        {0.0},
        // n=1
        {0.0, 0.0},
        // n=2
        {-0.484165371736E-03, 0.0, 0.243914352398E-06},
        // n=3
        {0.957161207093E-06, 0.202998882184E-05, 0.904627768605E-06, -0.721072657057E-06},
        // n=4
        {0.539873863789E-06, 0.350392272684E-06, -0.473670750654E-06, 0.662671572540E-06, 0.169994363185E-06}
    };
    
    private static final double[][] S_COEFFICIENTS = {
        // n=0
        {0.0},
        // n=1
        {0.0, 0.0},
        // n=2
        {0.0, 0.140016683654E-05, 0.0},
        // n=3
        {0.0, 0.248200415766E-06, -0.619005005177E-06, 0.0},
        // n=4
        {0.0, 0.662671572540E-06, 0.169994363185E-06, -0.473670750654E-06, 0.0}
    };
    
    /**
     * Calculate complete geoid height (using EGM2008 model)
     * 
     * @param latitude latitude in degrees
     * @param longitude longitude in degrees
     * @param maxDegree maximum degree (recommended not exceeding 360)
     * @return geoid height in meters
     */
    public static double calculateEGM2008Complete(double latitude, double longitude, int maxDegree) {
        double latRad = Math.toRadians(latitude);
        double lonRad = Math.toRadians(longitude);
        
        double geoidHeight = 0.0;
        
        // Calculate spherical harmonic expansion
        for (int n = 0; n <= maxDegree; n++) {
            for (int m = 0; m <= n; m++) {
                // Calculate Legendre function
                double Pnm = calculateLegendreFunction(n, m, latRad);
                
                // Calculate C coefficient terms
                if (m <= n && n < C_COEFFICIENTS.length && m < C_COEFFICIENTS[n].length) {
                    double Cnm = C_COEFFICIENTS[n][m];
                    if (m == 0) {
                        geoidHeight += Cnm * Pnm;
                    } else {
                        geoidHeight += Cnm * Pnm * Math.cos(m * lonRad);
                    }
                }
                
                // Calculate S coefficient terms
                if (m > 0 && m <= n && n < S_COEFFICIENTS.length && m < S_COEFFICIENTS[n].length) {
                    double Snm = S_COEFFICIENTS[n][m];
                    geoidHeight += Snm * Pnm * Math.sin(m * lonRad);
                }
            }
        }
        
        // Apply Earth's average radius and gravity constant
        double earthRadius = 6378137.0; // meters
        double gravityConstant = 3.986004418E14; // m³/s²
        
        return geoidHeight * earthRadius / gravityConstant;
    }
    
    /**
     * Calculate Legendre function P(n,m)
     * 
     * @param n degree
     * @param m order
     * @param theta colatitude (radians)
     * @return Legendre function value
     */
    private static double calculateLegendreFunction(int n, int m, double theta) {
        if (n < 0 || m < 0 || m > n) {
            return 0.0;
        }
        
        if (n == 0) {
            return 1.0;
        }
        
        if (n == 1) {
            if (m == 0) {
                return Math.cos(theta);
            } else if (m == 1) {
                return Math.sin(theta);
            }
        }
        
        if (n == 2) {
            if (m == 0) {
                return (3.0 * Math.pow(Math.cos(theta), 2) - 1.0) / 2.0;
            } else if (m == 1) {
                return 3.0 * Math.cos(theta) * Math.sin(theta);
            } else if (m == 2) {
                return 3.0 * Math.pow(Math.sin(theta), 2);
            }
        }
        
        // For higher degrees, use recursive formula
        return calculateLegendreRecursive(n, m, theta);
    }
    
    /**
     * Calculate Legendre function using recursive formula
     */
    private static double calculateLegendreRecursive(int n, int m, double theta) {
        double x = Math.cos(theta);
        double y = Math.sin(theta);
        
        // Initialize P(m,m) and P(m+1,m)
        double Pmm = 1.0;
        if (m > 0) {
            Pmm = (2 * m - 1) * y * Pmm;
        }
        
        if (n == m) {
            return Pmm;
        }
        
        double Pmp1m = (2 * m + 1) * x * Pmm;
        if (n == m + 1) {
            return Pmp1m;
        }
        
        // Recursively calculate P(n,m)
        double Pnm = 0.0;
        for (int k = m + 2; k <= n; k++) {
            Pnm = ((2 * k - 1) * x * Pmp1m - (k + m - 1) * Pmm) / (k - m);
            Pmm = Pmp1m;
            Pmp1m = Pnm;
        }
        
        return Pnm;
    }
    
    /**
     * Calculate geoid height gradient (for analyzing change rate)
     * 
     * @param latitude latitude in degrees
     * @param longitude longitude in degrees
     * @param delta small change amount for gradient calculation (degrees)
     * @return gradient vector [dN/dlat, dN/dlon]
     */
    public static double[] calculateGeoidGradient(double latitude, double longitude, double delta) {
        double N1 = calculateEGM2008Complete(latitude + delta, longitude, 20);
        double N2 = calculateEGM2008Complete(latitude - delta, longitude, 20);
        double N3 = calculateEGM2008Complete(latitude, longitude + delta, 20);
        double N4 = calculateEGM2008Complete(latitude, longitude - delta, 20);
        
        double dNdlat = (N1 - N2) / (2 * delta);
        double dNdlon = (N3 - N4) / (2 * delta);
        
        return new double[]{dNdlat, dNdlon};
    }
    
    /**
     * Calculate geoid height curvature (for analyzing local change characteristics)
     * 
     * @param latitude latitude in degrees
     * @param longitude longitude in degrees
     * @param delta small change amount for curvature calculation (degrees)
     * @return curvature matrix [[d²N/dlat², d²N/dlatdlon], [d²N/dlatdlon, d²N/dlon²]]
     */
    public static double[][] calculateGeoidCurvature(double latitude, double longitude, double delta) {
        double N00 = calculateEGM2008Complete(latitude, longitude, 20);
        double Npp = calculateEGM2008Complete(latitude + delta, longitude + delta, 20);
        double Npm = calculateEGM2008Complete(latitude + delta, longitude - delta, 20);
        double Nmp = calculateEGM2008Complete(latitude - delta, longitude + delta, 20);
        double Nmm = calculateEGM2008Complete(latitude - delta, longitude - delta, 20);
        
        double d2Ndlat2 = (Npp + Npm + Nmp + Nmm - 4 * N00) / (4 * delta * delta);
        double d2Ndlon2 = (Npp + Npm + Nmp + Nmm - 4 * N00) / (4 * delta * delta);
        double d2Ndlatdlon = (Npp - Npm - Nmp + Nmm) / (4 * delta * delta);
        
        return new double[][]{
            {d2Ndlat2, d2Ndlatdlon},
            {d2Ndlatdlon, d2Ndlon2}
        };
    }
    
    /**
     * Batch calculate geoid height for multiple coordinate points
     * 
     * @param coordinates coordinate point list
     * @param maxDegree maximum degree
     * @return geoid height result list
     */
    public static List<GeoidResult> calculateBatch(List<Coordinate> coordinates, int maxDegree) {
        List<GeoidResult> results = new ArrayList<>();
        
        for (Coordinate coord : coordinates) {
            try {
                double geoidHeight = calculateEGM2008Complete(coord.latitude, coord.longitude, maxDegree);
                double[] gradient = calculateGeoidGradient(coord.latitude, coord.longitude, 0.001);
                double[][] curvature = calculateGeoidCurvature(coord.latitude, coord.longitude, 0.001);
                
                GeoidResult result = new GeoidResult(coord, geoidHeight, gradient, curvature);
                results.add(result);
                
            } catch (Exception e) {
                System.err.println("Failed to calculate coordinate (" + coord.latitude + ", " + coord.longitude + "): " + e.getMessage());
            }
        }
        
        return results;
    }
    
    /**
     * Geoid height calculation result class
     */
    public static class GeoidResult {
        public final Coordinate coordinate;
        public final double geoidHeight;
        public final double[] gradient;
        public final double[][] curvature;
        
        public GeoidResult(Coordinate coordinate, double geoidHeight, double[] gradient, double[][] curvature) {
            this.coordinate = coordinate;
            this.geoidHeight = geoidHeight;
            this.gradient = gradient;
            this.curvature = curvature;
        }
        
        @Override
        public String toString() {
            return String.format("Coordinate(%.4f, %.4f): GeoidHeight=%.6fm, Gradient=[%.6f, %.6f]", 
                coordinate.latitude, coordinate.longitude, geoidHeight, gradient[0], gradient[1]);
        }
    }
    
    /**
     * Coordinate point class
     */
    public static class Coordinate {
        public final double latitude;
        public final double longitude;
        
        public Coordinate(double latitude, double longitude) {
            this.latitude = latitude;
            this.longitude = longitude;
        }
    }
}
